packages = c('rgdal', 'sf', 'tmap', 'tidyverse', 'sp', 'rgeos','maptools', 'raster', 'spatstat', 'tmaptools', 'spdep', 'OpenStreetMap', 'ggpubr', 'SpatialPosition', 'SpatialAcc', 'dplyr', 'shinycssloaders','plotly', 'shinythemes')
for (p in packages){
    if (!require(p, character.only = T)){
        install.packages(p)
    }
    library(p, character.only = T)
}


# Population Data
# columns = planning_area, subzone, elderly_count, total_count
popdata <- read_csv('data/planning-area-subzone-age-group-sex-and-type-of-dwelling-june-2011-2019.csv')
popdata2019 <- popdata %>%
  filter(year == 2019) %>%
  # filter(age_group == "65_to_69" | age_group == "70_to_74" | age_group == '75_to_79' | age_group == '80_to_84' | age_group == '85_to_89' | age_group == '90_and_over') %>%
  group_by(planning_area, subzone, age_group) %>%
  summarise(count = sum(resident_count)) %>%
  ungroup() %>%
  spread(age_group, count) %>%
  mutate(elderly_count = `65_to_69` + `70_to_74` + `75_to_79` + `80_to_84` + `85_to_89` + `90_and_over`) %>%
  mutate(total_count = rowSums(.[3:21])) %>%
  dplyr::select(planning_area, subzone, elderly_count, total_count)
popdata2019 <- mutate_at(popdata2019, .vars = c("subzone", "planning_area"), .funs=toupper)

# Eldercare Data
eldercare_sf <- st_read(dsn='data', layer='ELDERCARE')
## Reading layer `ELDERCARE' from data source `D:\IS415Team2\codes\data' using driver `ESRI Shapefile'
## Simple feature collection with 133 features and 18 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 14481.92 ymin: 28218.43 xmax: 41665.14 ymax: 46804.9
## projected CRS:  SVY21
eldercare_sf <- eldercare_sf %>%
    mutate(label = "Eldercare centres") 
eldercare_sf <- st_transform(eldercare_sf, 3414)
st_crs(eldercare_sf)
## Coordinate Reference System:
##   User input: EPSG:3414 
##   wkt:
## PROJCRS["SVY21 / Singapore TM",
##     BASEGEOGCRS["SVY21",
##         DATUM["SVY21",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4757]],
##     CONVERSION["Singapore Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",1.36666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",103.833333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",28001.642,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",38744.572,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Singapore"],
##         BBOX[1.13,103.59,1.47,104.07]],
##     ID["EPSG",3414]]
eldercare_sp <- as_Spatial(eldercare_sf)
eldercare_spatialpoint <- as(eldercare_sp, "SpatialPoints")

# Silver Infocomm Data
infocomm_sf <- st_read(dsn='data', layer='SILVERINFOCOMM')
## Reading layer `SILVERINFOCOMM' from data source `D:\IS415Team2\codes\data' using driver `ESRI Shapefile'
## Simple feature collection with 101 features and 18 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 12024.87 ymin: 28385.78 xmax: 41300.53 ymax: 47763.05
## projected CRS:  SVY21
infocomm_sf <- infocomm_sf %>%
    mutate(label = "Silver Infocomm Junc")
infocomm_sf <- st_transform(infocomm_sf, 3414)
st_crs(infocomm_sf)
## Coordinate Reference System:
##   User input: EPSG:3414 
##   wkt:
## PROJCRS["SVY21 / Singapore TM",
##     BASEGEOGCRS["SVY21",
##         DATUM["SVY21",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4757]],
##     CONVERSION["Singapore Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",1.36666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",103.833333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",28001.642,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",38744.572,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Singapore"],
##         BBOX[1.13,103.59,1.47,104.07]],
##     ID["EPSG",3414]]
infocomm_sp <- as_Spatial(infocomm_sf)
infocomm_spatialpoint <- as(infocomm_sp, "SpatialPoints")

# Chas Clinic Data
chas_sf <- st_read(dsn='data/chas-clinics-kml.kml')
## Reading layer `MOH_CHAS_CLINICS' from data source `D:\IS415Team2\codes\data\chas-clinics-kml.kml' using driver `KML'
## Simple feature collection with 1167 features and 2 fields
## geometry type:  POINT
## dimension:      XYZ
## bbox:           xmin: 103.5818 ymin: 1.016264 xmax: 103.9903 ymax: 1.456037
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84
chas_sf <- chas_sf %>%
    mutate(label = "Chas Clinics") %>%
    mutate(capacity = 1)
chas_sf <- st_transform(chas_sf, 3414)
st_crs(chas_sf)
## Coordinate Reference System:
##   User input: EPSG:3414 
##   wkt:
## PROJCRS["SVY21 / Singapore TM",
##     BASEGEOGCRS["SVY21",
##         DATUM["SVY21",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4757]],
##     CONVERSION["Singapore Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",1.36666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",103.833333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",28001.642,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",38744.572,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Singapore"],
##         BBOX[1.13,103.59,1.47,104.07]],
##     ID["EPSG",3414]]
chas_sp <- as_Spatial(chas_sf)
chas_spatialpoint <- as(chas_sp, "SpatialPoints")

# MP14 Subzone Data
mpsz_sf <- st_read(dsn='data', layer='MP14_SUBZONE_WEB_PL')
## Reading layer `MP14_SUBZONE_WEB_PL' from data source `D:\IS415Team2\codes\data' using driver `ESRI Shapefile'
## Simple feature collection with 323 features and 15 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
## projected CRS:  SVY21
mpsz_sf <- mpsz_sf %>%
    dplyr::select(SUBZONE_N, PLN_AREA_N, REGION_N, X_ADDR, Y_ADDR, SHAPE_Leng, SHAPE_Area, geometry)
mpsz_sf <- st_transform(mpsz_sf, 3414)
st_crs(mpsz_sf)
## Coordinate Reference System:
##   User input: EPSG:3414 
##   wkt:
## PROJCRS["SVY21 / Singapore TM",
##     BASEGEOGCRS["SVY21",
##         DATUM["SVY21",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4757]],
##     CONVERSION["Singapore Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",1.36666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",103.833333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",28001.642,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",38744.572,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Singapore"],
##         BBOX[1.13,103.59,1.47,104.07]],
##     ID["EPSG",3414]]
eldercare_sf <- st_join(eldercare_sf, mpsz_sf, join=st_intersects)
infocomm_sf <- st_join(infocomm_sf, mpsz_sf, join=st_intersects)
chas_sf <- st_join(chas_sf, mpsz_sf, join=st_intersects)

# HDB Data
hdb <- read_csv('data/hdb_data.csv')
hdb <- hdb %>% 
  mutate(total_count=rowSums(.[2:11])) %>%
  dplyr::select(Postcode, X, Y, Latitude, Longitude, total_count)
hdb$Postcode <- as.character(hdb$Postcode)
hdb_sf <- st_as_sf(hdb, coords=c('X', 'Y'), crs='EPSG:3414')

# SG Coastal Outline
sg <- readOGR(dsn = "data", layer="CostalOutline")
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\IS415Team2\codes\data", layer: "CostalOutline"
## with 60 features
## It has 4 fields
sg_spatialpoint <- as(sg, "SpatialPolygons")

# Compute elderly_density in mpsz
mpsz_sf <- left_join(mpsz_sf, popdata2019, by=c('PLN_AREA_N' = 'planning_area', 'SUBZONE_N' = 'subzone'))
mpsz_sf <- mpsz_sf %>%
  dplyr::select(SUBZONE_N, PLN_AREA_N, REGION_N, elderly_count, total_count, geometry) %>%
  mutate(elderly_proportion = elderly_count / total_count) %>%
  mutate_if(is.numeric, ~replace(., is.nan(.), 0))
mpsz_sp <- as_Spatial(mpsz_sf)
mpsz_spatialpoint <- as(mpsz_sp, "SpatialPolygons")

# Add count of facilities per subzone  into mpsz_demand
mpsz_demand <- mpsz_sf
mpsz_demand$`chas_count` <- lengths(st_intersects(mpsz_sf,chas_sf))
mpsz_demand$`eldercare_count` <- lengths(st_intersects(mpsz_sf,eldercare_sf))
mpsz_demand$`infocomm_count` <- lengths(st_intersects(mpsz_sf,infocomm_sf))

# Join HDB with MP14 Data
hdb_mpsz <- st_join(mpsz_sf, hdb_sf, join=st_intersects) %>%
  dplyr::select(SUBZONE_N, PLN_AREA_N, REGION_N, elderly_proportion, Postcode, total_count.y, geometry) %>%
  rename(resident_count = total_count.y) %>%
  mutate(elderly_count = resident_count * elderly_proportion) %>%
  mutate_if(is.numeric, ~replace(., is.na(.), 0))

#Finding area of each subzone
mpsz_demand$Area <- mpsz_demand %>%
    st_area()

#Calculating density of Chas, Eldercare, Infocomm and Elderly Pop
mpsz_demand <- mpsz_demand %>%
    mutate(`Elderly_Density` = (`elderly_count`/ Area) * 1000000) %>%
    mutate(`Chas_Density` = `elderly_count` / `chas_count`) %>%
    mutate(`Eldercare_Density` = `elderly_count` / `eldercare_count`) %>%
    mutate(`Infocomm_Density` = `elderly_count` / `infocomm_count`)

#Converting mpsz_demand to SP object
mpsz_demand_sp <- as_Spatial(mpsz_demand)
mpsz_demand_spatialpoint <- as(mpsz_demand_sp, "SpatialPolygons")

#reading OSM 
sg_osm <- read_osm(mpsz_spatialpoint, ext=1.3)

#Extracting unique planning areas
uniqPlanningAreas <- mpsz_demand[2]
uniqPlanningAreas <- st_set_geometry(uniqPlanningAreas, NULL)
uniqPlanningAreas <- unique(uniqPlanningAreas)

Quadrat Test

CHAS Clinics

owin <- as(sg, "owin")
spatialpoint <- as(sg, "SpatialPolygons")
chas_ppp <- as(chas_spatialpoint, "ppp")
chas_ppp_jit <- rjitter(chas_ppp, retry=TRUE, nsim=1, drop=TRUE)
ppp <- chas_ppp_jit[owin]  
ppp.km <- rescale(ppp, 1000, "km")
qt <- quadrat.test(ppp.km, nx=20, ny=15)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
qt
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp.km
## X2 = 1976.7, df = 184, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 185 tiles (irregular windows)
plot(ppp.km)
plot(qt, add=TRUE, cex=.1)

quadrat.test(ppp.km, nx=20, ny=15, method='M', nsim=2999)
## 
##  Conditional Monte Carlo test of CSR using quadrat counts
##  Test statistic: Pearson X2 statistic
## 
## data:  ppp.km
## X2 = 1976.7, p-value = 0.0006667
## alternative hypothesis: two.sided
## 
## Quadrats: 185 tiles (irregular windows)

Eldercare

owin <- as(sg, "owin")
spatialpoint <- as(sg, "SpatialPolygons")
eldercare_ppp <- as(eldercare_spatialpoint, "ppp")
eldercare_ppp_jit <- rjitter(eldercare_ppp, retry=TRUE, nsim=1, drop=TRUE)
ppp <- eldercare_ppp_jit[owin]  
ppp.km <- rescale(ppp, 1000, "km")
qt <- quadrat.test(ppp.km, nx=20, ny=15)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
qt
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp.km
## X2 = 495.78, df = 184, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 185 tiles (irregular windows)
plot(ppp.km)
plot(qt, add=TRUE, cex=.1)

quadrat.test(ppp.km, nx=20, ny=15, method='M', nsim=2999)
## 
##  Conditional Monte Carlo test of CSR using quadrat counts
##  Test statistic: Pearson X2 statistic
## 
## data:  ppp.km
## X2 = 495.78, p-value = 0.007333
## alternative hypothesis: two.sided
## 
## Quadrats: 185 tiles (irregular windows)

Silver Infocomm Junctions

owin <- as(sg, "owin")
spatialpoint <- as(sg, "SpatialPolygons")
infocomm_ppp <- as(infocomm_spatialpoint, "ppp")
infocomm_ppp_jit <- rjitter(infocomm_ppp, retry=TRUE, nsim=1, drop=TRUE)
ppp <- infocomm_ppp_jit[owin]  
ppp.km <- rescale(ppp, 1000, "km")
qt <- quadrat.test(ppp.km, nx=20, ny=15)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
qt
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  ppp.km
## X2 = 317.48, df = 184, p-value = 7.149e-09
## alternative hypothesis: two.sided
## 
## Quadrats: 185 tiles (irregular windows)
plot(ppp.km)
plot(qt, add=TRUE, cex=.1)

quadrat.test(ppp.km, nx=20, ny=15, method='M', nsim=2999)
## 
##  Conditional Monte Carlo test of CSR using quadrat counts
##  Test statistic: Pearson X2 statistic
## 
## data:  ppp.km
## X2 = 317.48, p-value = 0.05133
## alternative hypothesis: two.sided
## 
## Quadrats: 185 tiles (irregular windows)

Gi* Statistics - Hotspot & Coldspot Analysis

CHAS

#coverting to sp polygons
mpsz_demand_selected_sp <- as_Spatial(mpsz_demand)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum SVY21 in CRS definition
coords <- coordinates(mpsz_demand_selected_sp)
k1 <- knn2nb(knearneigh(coords))
k1dists <- unlist(nbdists(k1, coords, longlat=FALSE))
maxdist <- ceiling(max(k1dists))

#Calculating d fixed-distance weight matrix
dnb <- dnearneigh(coordinates(mpsz_demand_selected_sp), 0, maxdist, longlat = FALSE)

#plot(mpsz_demand_selected_sp, border = 'lightgrey')
#plot(dnb, coordinates(mpsz_demand_selected_sp), add=TRUE, col='red')

dnb_lw <- nb2listw(dnb, style = 'B')

fips <- order(mpsz_demand_selected_sp$SUBZONE_N)
gi.fixed <- localG(mpsz_demand_selected_sp$chas_count, dnb_lw)
gi.fixed
##   [1] -1.097895747 -1.750276041 -1.224729244 -1.649028543 -2.064212720
##   [6] -2.104086969 -1.901554875 -1.369558330 -2.593432159 -2.101292745
##  [11] -2.747198138 -1.546764134 -0.443062227 -1.503022100 -1.619944188
##  [16] -2.883406367 -1.032641641 -0.729049655 -1.034345429 -2.321395295
##  [21] -1.712554693 -1.380816216 -1.377488708 -2.354621229 -1.945165127
##  [26] -1.391393537 -1.480648481 -1.893688362 -2.471743794 -2.791361849
##  [31] -1.341978208 -1.900254240 -1.853839054 -1.577809366 -2.066549368
##  [36] -1.635487565 -1.350746913 -1.433679257 -1.406179580 -1.903140385
##  [41] -0.954469665 -1.845231592 -1.176629161  0.597843958 -0.951053930
##  [46] -0.578373599 -1.238172440 -1.201281267 -2.299024629  2.470425585
##  [51] -1.925531529 -1.102744323 -1.008256238 -0.222250411 -1.232462559
##  [56] -2.129614118 -0.956845194 -0.109051056 -0.061443232 -1.012196402
##  [61] -2.099018138 -1.401456246 -1.924088906 -1.804893203 -0.131215798
##  [66]  0.208445774 -1.785318927  0.607510237  1.244966241 -1.766319485
##  [71] -2.214459555 -2.142616072 -0.773761038 -1.198891758 -1.547837882
##  [76] -0.519272574  2.963033961 -1.518153429 -1.708842703 -0.884783701
##  [81] -0.524114540  1.206227662 -1.004389124 -0.460550974  0.170551313
##  [86] -1.448377090 -1.275613027 -1.453325491 -2.510166044 -0.775924961
##  [91] -0.781738019  3.867165618 -1.243749916 -1.364920993 -1.561403798
##  [96] -0.164861841 -0.400156963 -0.663477304  0.091373939 -0.181244996
## [101] -2.042811789 -1.096868369  2.263711543 -0.610344516  1.035179183
## [106] -0.461078963 -0.835681162  3.359510264  0.403639797 -0.092678780
## [111] -1.760848295 -1.149713099  1.929301293 -0.691916596 -0.689440003
## [116] -0.991327362  3.509434389  2.492025187 -1.783106419 -0.020486411
## [121] -0.316626207 -0.748973304 -0.640682356 -1.340390173 -1.759022989
## [126]  0.298219769 -0.487173476 -1.113412729 -1.707779143 -1.560725383
## [131]  0.150243991  1.006615778  0.444707511 -1.215477121 -0.700423409
## [136] -0.753898934  3.692363911  1.656606566  0.519779890 -1.767771086
## [141]  0.581051086  0.331250221 -0.744045789 -0.854269586  0.727213838
## [146] -0.752342166 -0.353957904 -0.135514849  4.816891117 -0.484375283
## [151] -0.665464144 -0.660591092  3.650728407 -0.543479673 -0.823613899
## [156] -0.804595171  2.226565860 -0.433703572  0.341962401 -0.605707770
## [161]  0.351531512  1.996550070  1.464655849 -1.770583299 -0.660028257
## [166] -0.615901112  2.701726336 -1.030122195 -1.585272657  2.726082529
## [171] -0.997068570 -0.378124714 -0.667194145 -0.656505414 -0.631230573
## [176]  2.116287814 -0.779633114 -1.191791490  4.438419343  1.797010574
## [181]  0.189959224  3.816613245 -0.553424696  4.109907249 -0.377496888
## [186]  0.410182962 -0.446047300 -0.951574046  2.830403238  3.762186860
## [191] -1.539150633 -0.407622050 -1.033057370  0.611124624 -1.118413741
## [196] -1.049475304  2.104382260  0.542541079  2.915818858 -0.024313786
## [201]  3.307348449  5.154034976  5.059951346  4.959442232  1.848243676
## [206] -0.485231877  0.018343884  4.603710056  3.854585776  0.180737056
## [211] -0.970553433 -0.304491459  3.617087201  3.115111856  0.055083410
## [216]  3.646815817  1.737154398  2.964538659  0.259977405  0.229945954
## [221]  2.627725713 -0.209747847 -0.094832404 -0.175113724  2.041998578
## [226]  3.288052444  2.060931592  2.009458242  2.675943006  2.038377353
## [231]  3.558856844  2.845958537 -0.777569553 -0.371987115  3.831941212
## [236]  0.262704797 -0.438097967  1.415084478  2.763624792  1.068759475
## [241]  1.004496532  1.579760687  4.983842000  2.442161253 -0.482342250
## [246]  2.785630534  3.927556462  1.912012549  1.940969274  1.686096189
## [251]  2.450373117 -0.270167897  3.059914933  3.734143532 -0.165299297
## [256]  4.341412597 -0.561380812  3.547369056  0.218680739  3.427092814
## [261]  1.260432448  1.425562372  3.149260431  0.455004649  3.398295195
## [266]  1.651349191  1.848801277 -0.016203023  3.866855809  2.143393269
## [271]  1.812849398  2.886512051  1.697836081  3.812651631  1.120558764
## [276]  2.175207367  0.286056637  0.068167348  0.915814740  1.666586860
## [281] -0.450310423 -0.370101867 -0.277393718  1.166187367  0.619002233
## [286]  1.351780797 -0.420218122 -0.168659987 -0.257462928 -1.045887948
## [291] -0.670982576 -0.049261540 -0.383978403 -0.402871897  0.163324809
## [296] -0.578440303 -0.517922395 -0.151407320 -0.402871897 -0.214008828
## [301] -0.130580152 -0.729049655 -0.496288425  1.284387953  0.121629718
## [306] -0.440804399 -0.169385777 -0.565079733  1.045040426  1.324366000
## [311]  0.276377584  1.482616983  1.164900538  1.141327426  0.306010504
## [316]  1.334154568  0.642349641 -0.292232586  1.421088188 -0.106483069
## [321] -0.454036964  0.009388915 -0.155004360
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = mpsz_demand_selected_sp$chas_count, listw = dnb_lw)
## attr(,"class")
## [1] "localG"
chas.gi <- cbind(mpsz_demand_selected_sp, as.matrix(gi.fixed))
names(chas.gi)[15] <- "gstat"
#colnames(chas.gi)[23] <- "gstat"

tm_shape(chas.gi) +
  tm_fill(col = "gstat", 
          style = "pretty",
          palette="-RdBu",
          title = "local Gi") +
  tm_borders(alpha = 0.5)
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## Variable(s) "gstat" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Silver Infocomm

#coverting to sp polygons
mpsz_demand_selected_sp <- as_Spatial(mpsz_demand)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum SVY21 in CRS definition
coords <- coordinates(mpsz_demand_selected_sp)
k1 <- knn2nb(knearneigh(coords))
k1dists <- unlist(nbdists(k1, coords, longlat=FALSE))
maxdist <- ceiling(max(k1dists))

#Calculating d fixed-distance weight matrix
dnb <- dnearneigh(coordinates(mpsz_demand_selected_sp), 0, maxdist, longlat = FALSE)

#plot(mpsz_demand_selected_sp, border = 'lightgrey')
#plot(dnb, coordinates(mpsz_demand_selected_sp), add=TRUE, col='red')

dnb_lw <- nb2listw(dnb, style = 'B')

fips <- order(mpsz_demand_selected_sp$SUBZONE_N)
gi.fixed <- localG(mpsz_demand_selected_sp$infocomm_count, dnb_lw)
gi.fixed
##   [1] -0.486612476 -1.376364058 -1.008135286 -1.205008596 -1.917588180
##   [6] -1.835112563 -1.515746348 -1.179188873 -1.760654330 -1.924128231
##  [11] -1.861020247 -1.234437804  0.149127773 -1.319506248 -1.319506248
##  [16] -1.645565831 -0.551036959 -0.389034578 -0.882050612 -1.641850457
##  [21] -1.064487068 -1.248137965 -1.467860028 -2.092487861 -1.771531599
##  [26] -1.275772612 -1.318327319 -1.518589515 -2.288893189 -1.882044464
##  [31] -1.290805083 -1.560752690 -1.738600475 -1.162327170 -1.936417346
##  [36] -1.403206798 -1.063604526 -0.977752383 -1.150709925 -1.261388659
##  [41] -0.934697775 -1.587150494 -1.063604526  0.421790648 -0.236677231
##  [46] -0.953129260 -1.252335640 -1.008135286 -1.928944488  2.088328803
##  [51] -1.374642034 -1.294917439 -1.294917439  0.003901357 -1.264393995
##  [56] -1.524727012 -1.294917439 -0.076361780 -0.039843614 -1.294917439
##  [61] -1.750221874 -1.190404776 -1.434030182 -1.487816736 -0.951978511
##  [66] -0.013825444 -1.502838754  0.343619428  0.709842797 -1.319506248
##  [71] -1.291053783 -1.768591752 -1.326050711 -0.114128677 -0.875378934
##  [76] -1.176441705  3.233741962 -0.821940607  0.210542781  0.282507712
##  [81] -1.261983618  0.239061681 -1.272460917 -1.133662674  0.063885256
##  [86]  0.217861017 -1.136513876 -0.994864430 -1.486036168 -1.443230893
##  [91]  1.153093188  3.547706765 -0.049266883  0.132227737 -1.319506248
##  [96] -0.179296261 -0.973171863 -1.198144937 -0.127106267 -0.004347985
## [101] -1.684449627 -1.283432859  1.846875306  1.247657145 -0.090438037
## [106]  1.153093188 -1.604659789  2.882880156 -0.294967230 -0.085358807
## [111] -0.494449761  0.085349263  0.721540310  0.783652921  0.687368524
## [116] -0.793907776  3.211302257  2.115279564 -1.777480801 -0.143145261
## [121] -1.176441705  0.503994319  0.245843706  0.017271964 -1.292241092
## [126]  0.115521520 -1.390330359  0.838833982 -0.844843424  0.826178498
## [131] -0.080845811  0.168626291 -0.110211408 -0.465143267  0.877980656
## [136] -0.342860900  2.567701077  0.582915024 -0.145219920  0.115445788
## [141] -0.136905521 -0.180199522  0.958002207  0.658991860 -0.387461980
## [146]  0.828357832  1.286702686 -1.401314417  3.481010487  1.286702686
## [151]  0.997662419  0.795022020  3.179633314  1.048620418 -1.500787116
## [156] -0.819098952  1.215230212  1.095332948 -0.456505456  1.001351757
## [161] -0.223477078  0.942163267  0.304532611 -1.773714843 -1.411987426
## [166]  0.750348145  2.951299975  0.635832044 -1.005233344  1.723438657
## [171]  0.769360026  0.854653912  1.057173571  1.044125256  0.747455913
## [176]  1.025748504  1.000442492 -0.132419064  2.921121419  0.894749185
## [181] -0.369523103  2.726888733  1.011143727  2.812836350  1.170642824
## [186] -0.521022855  1.420865341  0.167252189  2.620998241  2.328978335
## [191]  0.300711918  1.187411662  0.861530148  0.012105359 -1.647632551
## [196]  0.583073371  1.120556466 -0.223477078  3.284272272 -1.167513969
## [201]  2.246487378  3.135308815  3.208296650  3.094768036 -0.993580210
## [206]  0.601204966  1.222362162  2.786944939  2.848951994  0.271836140
## [211]  0.371328639 -0.312680521  3.727028771  3.552917404  1.515209891
## [216]  2.966972941 -0.574739535  1.900328919  1.347524120  1.501708413
## [221]  1.733797638  0.050968122  0.857044104  0.857044104 -0.528908873
## [226]  2.425195734 -0.589227639 -0.449999441  2.276903291 -0.729821079
## [231]  2.533257431  3.331309421  0.893907436 -0.875378934  2.824346690
## [236] -0.867118969  0.866905546 -0.607031413  1.862051396 -0.739070529
## [241] -1.162496176  0.336873331  2.992349884  1.564633673  0.569271437
## [246]  3.546558700  2.763942285 -0.183343531 -0.473621867 -0.782606833
## [251]  0.102466886  1.321905465  1.936338196  3.057574722  1.297014958
## [256]  2.924314950  0.085349263  3.455600764  1.673941810  1.951219237
## [261] -0.977351573 -1.067251071  2.180102726  1.234348476  2.763942285
## [266]  0.509209071  1.158429208  0.565629412  2.809967322 -0.821940607
## [271]  1.512947431  1.943021784 -0.422786180  2.563991079 -1.487632051
## [276] -0.348196319 -0.944508035 -1.874206604 -1.621224918 -0.778358810
## [281] -0.145078137 -1.046703609 -0.672850225 -1.303424879  2.320956289
## [286] -0.828211667 -1.696496133 -0.951635669  0.683429518 -1.635542723
## [291] -1.854800629 -1.420266926 -1.768359077 -1.480074459 -1.133527045
## [296] -1.529900947  0.300711918 -1.754732685 -1.480074459 -1.377921639
## [301] -1.788188205 -0.389034578 -1.650802421 -0.982254895 -1.558554170
## [306] -0.926921824 -0.281736893 -1.800177349 -1.132824855  0.209159055
## [311] -1.090483138  1.066992242  0.371267883 -0.807216324  0.915479515
## [316] -0.642394728 -1.614114182 -0.573319219  1.273907405 -1.318381795
## [321] -1.529900947 -1.325485892 -1.377921639
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = mpsz_demand_selected_sp$infocomm_count, listw = dnb_lw)
## attr(,"class")
## [1] "localG"
infocomm.gi <- cbind(mpsz_demand_selected_sp, as.matrix(gi.fixed))
names(infocomm.gi)[15] <- "gstat"
#colnames(infocomm.gi)[23] <- "gstat"

tm_shape(infocomm.gi) +
  tm_fill(col = "gstat", 
          style = "pretty",
          palette="-RdBu",
          title = "local Gi") +
  tm_borders(alpha = 0.5)
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## Variable(s) "gstat" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Eldercare

#coverting to sp polygons
mpsz_demand_selected_sp <- as_Spatial(mpsz_demand)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum SVY21 in CRS definition
coords <- coordinates(mpsz_demand_selected_sp)
k1 <- knn2nb(knearneigh(coords))
k1dists <- unlist(nbdists(k1, coords, longlat=FALSE))
maxdist <- ceiling(max(k1dists))

#Calculating d fixed-distance weight matrix
dnb <- dnearneigh(coordinates(mpsz_demand_selected_sp), 0, maxdist, longlat = FALSE)

#plot(mpsz_demand_selected_sp, border = 'lightgrey')
#plot(dnb, coordinates(mpsz_demand_selected_sp), add=TRUE, col='red')

dnb_lw <- nb2listw(dnb, style = 'B')

fips <- order(mpsz_demand_selected_sp$SUBZONE_N)
gi.fixed <- localG(mpsz_demand_selected_sp$eldercare_count, dnb_lw)
gi.fixed
##   [1]  2.093809425  1.318729340  1.698574954  0.614115327  0.674046348
##   [6]  0.531687060  1.339159707  1.612148831 -0.444692848 -0.234421411
##  [11] -0.971363307  1.115574781  2.067513040  1.372926746  1.372926746
##  [16] -1.847163555 -0.618544314 -0.436695075  1.236712454  0.350096505
##  [21]  0.653143550  1.917090489  2.067517888 -0.097971533  0.616302086
##  [26]  1.824586859  1.771542473  0.950192536 -0.236360393 -0.823440130
##  [31]  1.863194674  1.581660678  0.306931143  1.613853860  0.609955954
##  [36]  1.523766969  1.682282507  1.787074608  1.682282507  1.592337298
##  [41]  1.839819811  0.936735466  1.682282507  2.271290458  1.663511226
##  [46]  2.194658797  1.833906951  1.818781797 -0.004908287  1.987394780
##  [51]  0.559135889  2.139730569  2.139730569  2.547811491  1.868006495
##  [56]  0.244764736  2.139730569  2.190859221  2.429170865  2.139730569
##  [61]  0.699965243  1.931448968  1.130656597  0.356371372  2.241820468
##  [66]  2.463870810 -0.702763791  2.324121450  1.289110801  0.465214467
##  [71]  0.416195455  0.932081276  2.173898341  1.097832957 -0.982621316
##  [76]  2.141172103  2.238246073  0.707037179 -1.027893317 -0.877331911
##  [81]  2.159213134  2.615226252  2.309735762  2.191419676  2.529346174
##  [86] -1.023010008  1.663511226  1.363207076 -1.668089965  2.109325218
##  [91] -0.387470456  1.608045340 -0.921127080 -0.881287655  1.612148831
##  [96]  2.377478410  2.309843311  2.140707983  2.275338852  1.991886291
## [101]  0.584862277  2.224534580  1.801634120 -0.766656152  2.228337860
## [106] -0.387470456  2.116696707  1.048958800  2.301445676  2.174039173
## [111] -1.903839227 -0.789987361  2.563745958 -0.822752724 -0.737394625
## [116]  1.469655842  2.066337705  1.474232292  0.106271036  2.377652095
## [121]  2.377652095 -0.893708693 -0.735288835 -0.856209401 -0.560933211
## [126]  2.245021117  2.363732232 -0.469815627  0.289407435  0.462478944
## [131]  2.301445676  2.343107686  2.209240004  0.618884550 -1.023957135
## [136]  1.113501705  1.545293926  2.145310185  1.892011266 -1.828348667
## [141]  2.036154379  2.242716270 -0.797139415 -0.491527619  1.854734050
## [146] -1.355266687 -0.897644064  1.998960252  2.691078771 -0.897644064
## [151] -0.769768776 -0.424766489  1.533931833 -0.604199052  1.803324827
## [156]  1.243647025  1.545360967 -0.571478882  1.696651548 -1.023957135
## [161]  2.352169617  1.631669049  1.768207936  1.145722142  2.028878949
## [166] -0.677763415  1.940268440  0.523336583  0.564511093  1.367536388
## [171] -0.692841005 -0.845748361 -0.971187610 -1.113864123 -1.138058248
## [176]  0.932081276 -1.095495196 -1.668089965  2.818601813  1.032592479
## [181]  1.852829350  0.734585952 -0.985950335  2.662697137 -0.864375684
## [186]  0.963594686 -0.581080557 -0.231627044 -0.421498973  2.055646099
## [191] -0.952934588 -1.180051357 -0.735288835  0.851317257  1.714492812
## [196]  0.076187155  1.425944264  2.239237880  0.864663507  2.072479360
## [201]  0.358418356  2.797607307  2.806962143  2.723733544  0.289407435
## [206]  0.071683671 -0.509826981  3.007235385  1.802136107  0.742773366
## [211] -0.874602082  1.613853860  0.717103125  0.560179418 -0.332568742
## [216]  0.266997709  0.693068610  2.083517689  0.089514422  0.378612796
## [221]  0.272253000  0.399168079  0.270875752  0.270875752  1.057629860
## [226]  2.374192653  1.302198555  1.083785797  0.325579368 -0.014506531
## [231]  2.532334922  0.316796962 -0.832014308 -0.982621316  1.755142250
## [236]  0.784945513 -0.877331911  0.968921740  2.025745544  0.776694173
## [241]  0.319352005  1.391222669  2.812491699 -0.079194981  0.793992923
## [246]  0.160707278  0.856195800  1.028701092  0.974676599  0.718697669
## [251]  1.318808130 -0.836785982  1.926158338  0.311743554 -0.877331911
## [256]  1.823077978 -0.605523308  0.484029493  0.451002857  1.216328192
## [261]  0.141117129  0.594160932 -0.300799512  0.218927060 -0.056198180
## [266] -0.498120501 -0.664460833 -0.205382853  1.169024087  0.671026347
## [271] -0.398818264  1.355470314  0.764150938  0.992377247 -0.038016008
## [276] -0.110930387 -1.460641147 -0.208263843  0.207574817 -0.836785982
## [281] -1.235140358 -1.508093443 -0.918081219 -0.759286885  0.757378365
## [286] -0.892969903 -0.424766489 -1.127485479 -0.239381866 -1.513017996
## [291] -0.934132873 -1.047357068 -0.886146303 -1.038550292 -1.460641147
## [296] -1.101432806 -0.253440050 -1.127485479 -1.038550292 -0.909056676
## [301] -1.163171933 -0.436695075 -1.064905886 -0.360329507 -0.556049735
## [306] -1.211192661 -1.627693317 -0.858054530 -0.686613784 -0.571478882
## [311] -0.527921566 -0.715914034  0.746107306 -0.723876456  0.387238121
## [316] -0.723876456 -0.270815112 -1.064905886 -0.694276487 -1.462044259
## [321] -1.101432806 -0.842279246 -0.909056676
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = mpsz_demand_selected_sp$eldercare_count, listw = dnb_lw)
## attr(,"class")
## [1] "localG"
eldercare.gi <- cbind(mpsz_demand_selected_sp, as.matrix(gi.fixed))
names(eldercare.gi)[15] <- "gstat"
#colnames(eldercare.gi)[23] <- "gstat"

tm_shape(eldercare.gi) +
  tm_fill(col = "gstat", 
          style = "pretty",
          palette="-RdBu",
          title = "local Gi") +
  tm_borders(alpha = 0.5)
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## Variable(s) "gstat" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

Kernel Density Maps

CHAS Clinics

bindingbox <- st_bbox(sg)
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
owin <- as(sg, "owin")
spatialpoint <- as(sg, "SpatialPolygons")
chas_ppp <- as(chas_spatialpoint, "ppp")
chas_ppp_jit <- rjitter(chas_ppp, retry=TRUE, nsim=1, drop=TRUE)

ppp_chas <- chas_ppp_jit[owin]

selected_osm <- read_osm(bindingbox, ext=1.1)

msgtext <- "KDE Map"

ppp.km <- rescale(ppp_chas, 1000, "km")

kde <- density(ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
## Warning: point-in-polygon test had difficulty with 215 points (total score not 0
## or 1)
gridded_kde_bw <- as.SpatialGridDataFrame.im(kde)

kde_bw_raster <- raster(gridded_kde_bw)

projection(kde_bw_raster) <- crs("+init=EPSG:3414 +datum=WGS84 +units=km")

# Plot kernel density map on openstreetmap
tm_shape(selected_osm)+ 
  tm_rgb()+
  tm_shape(kde_bw_raster) + 
  tm_raster("v", alpha=0.5,  
            palette = "YlOrRd")

Silver Infocomm Junctions

bindingbox <- st_bbox(sg)
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
owin <- as(sg, "owin")
spatialpoint <- as(sg, "SpatialPolygons")
infocomm_ppp <- as(infocomm_spatialpoint, "ppp")
infocomm_ppp_jit <- rjitter(infocomm_ppp, retry=TRUE, nsim=1, drop=TRUE)

ppp_infocomm <- infocomm_ppp_jit[owin]

selected_osm <- read_osm(bindingbox, ext=1.1)

msgtext <- "KDE Map"

ppp.km <- rescale(ppp_infocomm, 1000, "km")

kde <- density(ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
## Warning: point-in-polygon test had difficulty with 215 points (total score not 0
## or 1)
## Warning: point-in-polygon test had difficulty with 215 points (total score not 0
## or 1)
gridded_kde_bw <- as.SpatialGridDataFrame.im(kde)

kde_bw_raster <- raster(gridded_kde_bw)

projection(kde_bw_raster) <- crs("+init=EPSG:3414 +datum=WGS84 +units=km")

# Plot kernel density map on openstreetmap
tm_shape(selected_osm)+ 
  tm_rgb()+
  tm_shape(kde_bw_raster) + 
  tm_raster("v", alpha=0.5,  
            palette = "YlOrRd")

Eldercare Services

bindingbox <- st_bbox(sg)
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
owin <- as(sg, "owin")
spatialpoint <- as(sg, "SpatialPolygons")
eldercare_ppp <- as(eldercare_spatialpoint, "ppp")
eldercare_ppp_jit <- rjitter(eldercare_ppp, retry=TRUE, nsim=1, drop=TRUE)

ppp_eldercare <- eldercare_ppp_jit[owin]

selected_osm <- read_osm(bindingbox, ext=1.1)

msgtext <- "KDE Map"

ppp.km <- rescale(ppp_eldercare, 1000, "km")

kde <- density(ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
## Warning: point-in-polygon test had difficulty with 215 points (total score not 0
## or 1)
## Warning: point-in-polygon test had difficulty with 215 points (total score not 0
## or 1)
gridded_kde_bw <- as.SpatialGridDataFrame.im(kde)

kde_bw_raster <- raster(gridded_kde_bw)

projection(kde_bw_raster) <- crs("+init=EPSG:3414 +datum=WGS84 +units=km")

# Plot kernel density map on openstreetmap
tm_shape(selected_osm)+ 
  tm_rgb()+
  tm_shape(kde_bw_raster) + 
  tm_raster("v", alpha=0.5,  
            palette = "YlOrRd")